/*
 * RumpfCO2Model.cpp
 *
 *  Created on: 26 Oct 2011
 *      Author: allan
 */

#include "RumpfCO2Model.h"

// GeoReact includes
#include "Geochemistry/ElectrolyteSolution.h"
#include "Utils/Algorithms.h"

const double RumpfCO2ActivityCoefficient(const ElectrolyteSolution& sol, double T, double P, const VectorXd& n)
{
	// The temperature in degrees Kelvin
	const double Tk = T + 273.15;
	
	// The ions in the electrolyte solution
	const vector<string>& ions = sol.GetIons();
	
	// The indexes of the ions Cl[-], Na[+], K[+], Ca[2+] and Mg[2+]
	const unsigned iCl = IndexOf( "Cl[-]", ions);
	const unsigned iNa = IndexOf( "Na[+]", ions);
	const unsigned iK  = IndexOf(  "K[+]", ions);
	const unsigned iCa = IndexOf("Ca[2+]", ions);
	const unsigned iMg = IndexOf("Mg[2+]", ions);
	
	// The stoichiometric molalities of the ions
	const VectorXd mi = sol.StoichiometricIonicMolalities(n);
	
	// The number of ions in the electrolyte solution
	const unsigned Ni = ions.size();
	
	// The stoichiometric molalities of the ions Cl[-], Na[+], K[+], Ca[2+] and Mg[2+]
	const double mCl = (iCl < Ni) ? mi[iCl] : 0.0;
	const double mNa = (iNa < Ni) ? mi[iNa] : 0.0;
	const double mK  = (iK  < Ni) ? mi[iK ] : 0.0;
	const double mCa = (iCa < Ni) ? mi[iCa] : 0.0;
	const double mMg = (iMg < Ni) ? mi[iMg] : 0.0;
	
	// Auxiliary temperature variables
	const double Tk1 = 1.0/Tk;
	const double Tk2 = Tk1 * Tk1;
	const double Tk3 = Tk2 * Tk1;
	
	// Pitzer's paramenters
	const double B0    = 0.254 - 76.82 * Tk1 - 10656.0 * Tk2 + 6312.0E3 * Tk3;
	const double Gamma = -0.0028;
	
	// The activity coefficient for aqueous CO2 computed with the model of Rumpf et al. (1994) extended for other chloride solutions by Spycher and Pruess (2005), Appendix A.2
	const double gCO2aq = std::exp(2 * B0 * (mNa + mK + 2*(mCa + mMg)) + 3 * Gamma * mCl * (mNa + mK + mCa + mMg));
	
	return gCO2aq;
}
